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Abstract. In this article, we present molecular dynamics study of the velocity autocorrelation function 
(VACF) of a Brownian particle. We compare the results of the simulation with the exact analytic predictions 
for a compressible fluid from [6] and an approximate result combining the predictions from hydrodynamics 
at short and long times. The physical quantities which determine the decay were determined from separate 
bulk simulations of the Lennard- Jones fluid at the same thermodynamic state point. We observe that the 
long-time regime of the VACF compares well the predictions from the macroscopic hydrodynamics, but 
the intermediate decay is sensitive to the viscoelastic nature of the solvent. 
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1 Introduction 

Brownian motion is the random motion of a particle, which 
is large compared to the solvent molecules, but is not 
of macroscopic size. It has become a paradigm in vari- 
ous branches of science and remains an active area of re- 
search among theoreticians and experimentalists. It is not 
only a preferred tool of theoretical modeling, but is also 
extensively used to probe microscopic environments in 
experiments [MPmir. A considerable effort is also spend 
in investigating transport properties of colloidal suspen- 
sions and complex fluids, primarily due to the relevance 
of such systems in industry, using molecular simulations 
of Brownian motion [22 13 18 14 8J. Such simulations al- 
ways invariably involve using discrete particles to explain 
continuum predictions of theory, and therefore, requires 
a clear understanding of the molecular and continuum 
regimes. 

In this article, using molecular dynamics simulation, 
we investigate the velocity autocorrelation function (VACF) 
of a Brownian particle. We choose a large system size so 
that the effect of finite-size of the simulation box is small. 
The internal degrees of freedom of the nanoparticle is also 
resolved in the simulation, such that stick boundary condi- 
tions apply on the surface of the particle [TS]. The erratic 
motion of a Brownian particle exhibits a far more rich 
behavior than predicted by the simple Langevin picture. 
In the continuum description, a fiuid is well described by 
the Stokes equation, with the particle dynamics coupled 
to the solvent through the imposed boundary conditions. 
The inadequacy of the Langevin picture in describing such 
erratic motion can be immediately seen from the VACF of 
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the Brownian particle. While a simple exponential decay 
is predicted at all times in the Langevin model, in real- 
ity the decay exhibits distinct features of both continuum, 
as well as discrete nature of the solvent. Accordingly, we 
classify the decay in three separate regimes, a short-time 
regimes - where molecular nature of the solvent plays a 
crucial role, an intermediate regime - governed by the in- 
terplay between sound propagation, vorticity diffusion and 
the viscoelasticity of the solvent, and a long-time regime 
where the VACF decays as a power law t^"^!"^ (d is the 
dimension of space) due to the development of the slow 
viscous patterns in the solvent |1|2I12|6|5] . We compare 
the results from the simulations with the exact predic- 
tions from hydrodynamics, and observe, that, while the 
decay of the VACF at long times compares well with the 
predictions from hydrodynamics, the intermediate decay 
is sensitive to the viscoelastic nature of the fluid and can 
not be explained by only considering the compressible na- 
ture of the fluid. 

Using a molecular dynamics simulation to investigate 
short and long-time dynamics of isothermal Brownian mo- 
tion is a non-trivial task. The key points in such simula- 
tions, are the identiflcation of relevant length and time 
scales. The two important length scales in the system are 
the simulation box size L and the radius of the particle 
R. In a typical molecular dynamics simulation periodic 
boundary conditions are imposed, implying that the dy- 
namics of the particle is effected by its periodic images. 
Since the strength of such finite size effect is determined 
by the ratio of the two length scales, R/L, we have the 
choice of a small R or large L. Both of these choices are 
unfortunately restricted. While the choice of L is solely 
determined by the computational resources at hand, the 
choice of R is determined by a number of factors. Ideally, 
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one would prefer a clear separation of the time scales in 
the simulation, in particular, the sonic time = i?/c and 
the vorticity diffusion time = fv^ both of which de- 
termine the decay of the VACF. A small choice of R does 
not resolve these time scales properly. 

The lower bound for R is determined by the Knudsen 
number, defined as the ratio of the mean free path to 
the characteristic length scale of the flow - typically, the 
diameter of the particle. The Knudsen number decides 
whether a continuum or a statistical mechanics description 
of the system is appropriate. Besides the length scales, the 
time scales involved in a Brownian motion range from the 
order of 10"^^ s to seconds. The simulation must also be 
able to resolve the various time scales in the problem, the 
smallest of which is the collision time of solvent molecules 
and the largest time scale is the colloid diffusion time, over 
which the colloid diffuses over its own radius. 

The remainder of the article is organized as follows. In 
Section [2] we explain our molecular dynamics simulation 
in brief. The comparison of the simulation results with the 
exact prediction from the theory is done in Section [Sj and 
an approximate result is presented in Section |4] on page[5j 



2 Molecular Dynamics Simulation 

In this section, we provide the details of our molecular 
dynamics simulation. To begin with, the natural choice of 
units is the Lennard- Jones reduced units, where length, 
time and energy in units of cr, t = \Jma'^ je and e. Through- 
out the article, the numerical values of the physical quanti- 
ties are given in reduced units, unless otherwise explicitly 
stated. 

Our model system is made of a simple Brownian par- 
ticle, with internal degress of freedom resolved, immersed 
in a Lennard- Jones solvent. The particles in the system 
interact via the Lennard-Jones interaction. 



U{r) = 4e 



(1) 



Additionally, the nearest-neighbor interaction of the atoms 
in the spherical cluster of the nanoparticle is the FENE 
interaction. 
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UpENEir) = --hRq log 



(2) 



where k = 30e/a^ is the spring constant and Rq = 1.5cr. In 
order to keep the finite-size effects from the image particles 
to a reasonable value, we choose a large system size with 
a total of 256000 particles. The initial configuration of the 
system was chosen to be a perfect FCC lattice with veloc- 
ities drawn from a Boltzmann distribution. The nanopar- 
ticle was obtained from a spherical cut of the FCC lattice. 
For the smallest size of the nanoparticle (i? = 3) there 
were 177 atoms in the cluster while for largest size of the 
nanoparticle (R = 5) we have 767 atoms in the cluster. 



The ratio L/R, which quantifies the finite-size in the sys- 
tem are 22.5795 and 13.5409 for the i? = 3 and i? = 5, 
respectively. 

The system was first equilibrated under NPT ensem- 
ble, with the system coupled to a thermostat and barostat, 
at a thermodynamic pressure of Pq = 0.01 and tempera- 
ture Tg = 0.75. For the implementation of the NPT en- 
semble we chose the Nose-Hoover equations of motion as 
modified by Melchionna|16j. 

Ti = Vi + a(fi - Ro), Pi = F,; - (a + 7)(f,; - Rq) 
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V = Wa. 



(3) 



The barostating variable a can be eliminated between 
the last two equations in Eq.(|3]), and the resulting equa- 
tions are then numerically integrated using Leap-Frog in- 
tegration scheme [21'. 

The molecular dynamics simulations were implemented 
on Graphics Processing Units (GPU) and is similar to 
those of Anderson et. al. [3] and Colberg |7j, more closely 
resembling the later in the construction of the Verlet list. 
We briefly describe our implementation in the following 
lines. 

We use the atom decomposition method in the simula- 
tion, for efficient parallel implementation of our MD code. 
Every particle in the simulation is assigned a thread, which 
is responsible for updating the coordinates and momenta 
of the particle. A GPU optimized cell list algorithm is used 
for construction of the Verlet list. For this purpose, the 
simulation domain is divided into cubes of size Tc, where 
Tc is the cutoff length scale for the Lennard-Jones poten- 
tial = 2.5a in the present application). The particles 
were first sorted into their respective cells using the paral- 
lel radixsort algorithm [7]. To construct of the Verlet list, 
the entries of 26 + 1 cells are copied to the shared memory. 
Every cell has an upper limit for the maximum number of 
entries, determined by the size of the shared memory on 
the GPU. This limitation also prevented us from simul- 
taneously copying the particle coordinates to the shared 
memory. For a given particle in a cell, an iterative search is 
made of the neighboring cells and the particle coordinates 
are read from a texture array. The stability and numerical 
accuracy of our molecular dynamics code was verified by 
outputting the total energy and total momentum of the 
system. With an integration time-step 6t = 0.001, trajecto- 
ries of 2 X 10'' steps (corresponding to a physical duration 
of 20 ns) were simulated for the data points. 



3 Velocity Autocorrelation Function 

The long time tails in the VACF can be explained by the 
generalized Langevin equation. 



(4) 
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together with the time dependent friction coefficient C(Oi 
where P is the momentum of the Brownian particle and 
^ is the random force acting on it. Multiplying Eq.Q by 
P(0), and taking an ensemble average, the velocity auto- 
correlation function for the Brownian particle becomes, 



Cit) = j, [\{t-t')C{t'), 
M Jo 

which, in the frequency space is written as, 

C(^) 



C{lu) = C(0) 



\ M J 



C{0)Y{u). 



(5) 



(6) 



The zero- frequency limit of Eq.(|6| produces the Stokes- 
Einstein relation D = k^TjC^Q. Transforming back to real 
time, the normalized VACF of the Brownian particle is 
given by 



C{t) 
CiO) 



r 



2^ 



Re[Y{u)]cos{ujt) 



(7) 



For an incompressible fluid, the frequency dependent 
friction coefficient is given by |12)24j : 

C(w) = 67r77oi?(l + R{-iuj/iyy/^ - {iujR^ /9iy)), (8) 

where, 770 and v are the steady-state dynamic and kine- 
matic viscosities of the solvent, respectively. The square- 
root singularity in Eq. ([s]) gives rise to the power law decay 
of the VACF at long times. Because of the incompressibil- 
ity condition, the equal-time value of the VACF suffers 
a discontinuity from the equipartition value of ksT/M 
to kBT/M*, with the effective mass M* given by M* = 
M + Mf/2. Mf is the mass of the displaced fluid. 

In simulations, however, the discontinuity is not ob- 
served due to the finite compressibility of the fluid [241416] . 
In a compressible solvent, the sound propagation occurs 
with a finite speed, and the solvent surrounding the col- 
loid is not instantly set in motion. A fraction of the en- 
ergy of the Brownian particle is thus spent in creating 
these sound waves. To account for the compressibility ef- 
fect of the solvent, we need to consider the Boussinesq 
force for unsteady motion in a compressible solvent ^24^ . 
For this purpose, we consider the frequency dependent 
friction coefficient presented in the works of Chow et. al. 
[6] . For simplicity, we follow the notations in [6 . In terms 
of the vorticity diffusion time t^, = B? jv and the sonic time 
Tc = i?/c, C,{uj) is written as 



X = i\/iar7v, y = iuJTc 
4 



i^T^ (^+(4/3)7?) 



-1/2 



CH = --TTiqRx'lil - y)Q + 2{x - 1)P], 



(9) 



where P and Q are functions of the dimcnsionless variables 
X and y. We refer the readers to 6 , for an explicit ex- 
pression of these functions. Substituting Eq.(|9]) in Eq.([6]), 
and subsequently using we obtain the normalized 

VACF in real time. 



The physical parameters which enter Eq.([9]) were de- 
termined from separate molecular dynamics simulations of 
the bulk Lennard- Jones fluid at the same thermodynamic 
state point. The shear viscosity rj and the bulk viscosity 
fi were estimated from the off-diagonal and diagonal com- 
ponents of the stress tensor and using the Green-Kubo 
formula. 



krJ Jo 



(10) 



^(*) = r^ rmt')5pmdt' (11) 



The infintie time limit of Eq.( 10 ) and Eq.( 11 ) provided the 



steady state values of the shear and bulf viscosity, ryo and 
fiQ, respectively. The adiabatic sound speed in the solvent 
was estimated from the relation 



7 

prnxr' 



(12) 



where 7 is the ratio of the specific heats Cp/Cy, p is the 
density of the solvent, m is the mass of the solvent particles 
and xt is the isothermal thermal compressibility. 

Before we compare the results of the molecular dynam- 
ics simulations with the predictions from hydrodynamics, 
the assumptions of the macroscopic theory needs to val- 
idated. The crucial assumption which enters the theory 
is the boundary condition on the surface of the particle. 
While Eq.([9| assumes a stick boundary condition on the 
particle surface, such an assumption may break down in a 
microscopic scale. To validate this we measured the fric- 
tion coefficient for different radius of a rough Brownian 
particle. In figure [l] we show this dependence of the fric- 
tion coefficient on the radius of the Brownian particle and 
compare it with predictions from Stokes law with stick 
and slip boundary condition. The reasonable agreement 
of the steady-state friction coefffcient with the Stokes law 
Co = GTrrjoR indicates that hydrodynamic boundary con- 
ditions applicable on the surface of the sphere are those 
of stick boundary conditions!^ Additionally, the Kundsen 
number (Kn) for the system was determined by measuring 
the mean-free path of the solvent from the ballistic regime 
of the solvent mean-square displacement. The Knudsen 
number determines whether a statistical mechanics or a 
continuum description is more appropriate for the sys- 
tem and for large values of Kn (typically Kn > 0.1), de- 
viations from the continuum description become relevant. 
The measured values of Kn were Kn = 0.02261 and Kn = 
0.01366, for i? = 3 and i? = 5, respectively, which indicates 
that the assumptions in the macroscopic hydrodynamics 
remains valid in the present scenario. 



^ When the radius of the Brownian particle is comparable 
to the size of the solvent particles, the Stokes-Einstein relation 
with standard stick or slip boundary conditions can break down 
and a non-standard boundary condition may be required [15) . 
However, for the particle sizes for which the velocity autocor- 
relation function has been investigated in the present article, 
stick boundary conditions are valid as depicted in Figure [l] 
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Fig. 1. The dependence of the friction coefficient for dif- 
ferent radius of the nanoparticle. The solid line is plot of 
Stokes law for stick boundary condition, Co = 67r?7oi?, while 
the dashed line shows the plot of Co = 47r?7oi?. The radius 



of the nanoparticle were estimated from the Eq.( 14) 
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Fig. 2. Velocity autocorrelation function of a Brownian 
particle with radius i? = 5 (□,■) and R = 3 (°, •) for 
two different sizes of the simulation box, L « 68 (□,o) and 
L« 85 (□,•). 



Moreover, due to the periodic boundary condition im- 
posed in the simulations, the data suffers from finite-size 
effects from the long-ranged flow field of the image parti- 
cles. We chose a sufficiently large simulation box, so that 
the artifact of the image particles are small. To substanti- 
ate this, we performed simualtions with two different box 
lengths, L Rj 68 and L « 87, the result of which is depicted 
in the Figure |2] The measured velocity autocorrelation 
functions does not exhibit pronounced finite-size effects, 
particulary in the intermediate regime of interest. 

To compare the reuslts from the simulation with the 
theoretical predictions for the VACF, the numerical eval- 
uation of Eq. ([? ) was first carried out with the steady state 
values of the shear viscosity. However, we observed that 
the intermediate decay of the VACF can not be accurately 
described by only treating the solvent as compressible and 
it was essential to consider the viscoelastic nature of the 
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Fig. 3. Decay of 1 - ?7(i)/?7o with time. Two distinct 
decay-time scales are seen. At short times, comparable 
with molecular collision time, the decay is faster compared 
to the intermediate regime where power-law decay of the 
stress autocorrelation function becomes important. In the 
time-scales of interest, the decay can be modeled as sim- 
ple exponential with decay time ti and T2. The two time 
constants ri (corresponding to the dashed line) and T2 
(corresponding to the dot-dashed line) obtained from the 
above plot are 0.114 and 0.390, respectively. The inset 
shows the time dependent viscosity normalized by ryp. 

solvent [23] . As already pointed out in [TT] , the interaction 
of a colloid in a viscoelastic solvent can be visualized by 
considering the colloid connected to the fluid by a dash- 
pot and spring in series. At times larger than t^, the vis- 
cous dissipation is represented by the dash-pot, while at 
shorter times the colloid interacts with the fluid via elas- 
tic forces. This caging effect is more pronounced when the 
ratio of the mass of the Brownian particle to the solvent 
particles is small [T7]. To this end, we model the Lennard- 
Jones solvent as Maxwell fluid with the frequency depen- 
dent viscosity given by 



fi{uj) = ?7o/(l -iwr). 



(13) 



where 770 is the steady-state shear viscosity of the solvent 
(lu = component). The relaxation time r is related to 
the infinite frequency shear modulus Goo as t = rjo jG 00 ■ 
For a Lennard- Jones solvent, the single exponential relax- 
ation in Eq.(13l ignores the algebraic decay of the stress 



autocorrelation at long times. This is also evident from 
the time dependent viscosity obtained from the simula- 
tion using Eq.([To|). To illustrate this more clearly, we plot 
the variation of the quantity 1 - ?7(i)/?7o with time in Fig. 
|3] The two distinct decays shown in Fig. |3] can be mod- 
elled as simple exponential with decay times ri and T2- 
To determine the VACF in the intermediate regime using 
Eq.0 and Eq.([9]) , we use T2 as the relaxation time for the 
model fluid. Since the colloid is made of discrete number 
of particles, the radius of the sphere R was determined 
from the radius of gyration using the relation 
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Fig. 4. Decay of the velocity autocorrelation function of 
a Brownian particle for two particle radius R = 5 (□) and 
i? = 3 (o). The solid-line is numerical evaluation of the 
normalized VACF using Eq.([7]) and Eq.Q together with a 
frequency dependent viscosity (Eq.( 13 )). The decay in the 
intermediate regime is sensitive to tEe viscoelastic nature 
of the solvent, and can not be accurately described using 
only the assumption of compressible fluid (see inset). The 
dashed line in the inset is the numerical evaluation with 
a constant viscosity. 



Physical 
quantities 
V 
M 
c 



Value used in nu- 
merical evaluation 

2.79961 

0.78561 

5.15444 



From bulk 
simulation 
2.79961 ±0.12963 
0.78556 ± 0.07442 
5.18256 ±0.05225 

Table 1. Comparison of the numerical values of the phys- 
ical quantities 77, /z and c, obtained from the bulk simu- 
lations of Lennard- Jones fluid and from best fit of the 
simualted VACF using Eq.(7]). We note, that the values of 
these parameters agree witmn statistical error bars. 



Fig. 5. Comparison of the theoretical prediction of VACF 
with simulation for two values of the radius of the parti- 
cle. The dot-dashed lines is the plot of addition of Eq.p^ 
and Eq.(16|. Even though the fit to the simulation data 



is poor in the short and intermediate regime, it still re- 
produces the qualitative nature of the decay. The physical 
parameters used in the plot are the same as presented in 
Table ID 



viscous patterns [9 , with usually an order of magnitude 
larger compared to t^. Assuming this time scale separa- 
tion, the complete decay of the VACF can be constructed 
by a simple addition of the VACF at short and long-time 
regime [TS]. 

In the long-time regime, following [I2[IS], the normal- 
ized VACF of a Browninan particle can be written as: 



M{V{t)V{0))/kBT: 



2a 

3pf 37r 



-r 



dx 



In the numerical evaluation of the VACF, we observed that 

the radius of the particle which gives a more accurate fit . 

to the data was close to i? + a/2, where a is the diameter T:{l-^crf - Aa2)/y 2{ai 
of the Lennard- Jones fiuid particles. The values of R used 
were 5.35 and 3.33 compared to the value of R + a/2 = 5.37 
and R+a/2 = 3.42, respectively. In table[l] we compare the 
numerical values of the physical quantities 77, fi and c for 
a bulk Lennard- Jones fluid and those which provided the 
best fit to the simulated velocity autocorrelation funtion 
using Eq. ([t]) . Additionally, the mass of the colloidal par- 
ticle was always taken as the reduced mass of the system. 



1 + aix + a2x'^ 

(15) 

(l/9)(7 - Ap,/pf) and = (1/81)(1 + 2p,/pfY 
In Eq.(|15), Pc is the density of the colloid and P/ is the 
density oTthe fluid. At i = 0, the integral gives the value 



with ai 



\l a\- 4(72 ) , and the right-hand 
side of Eq.(15), after simplification, becomes {2pc/ Pf)/{1+ 
2pc/pf), producing the well known discontinuity. The dis- 
continuity is quite easily removed when we take into ac- 
count the finite compressibility of the fiuid. To this end, 
we consider the general expression of Zwanzig |25j , 



M{V{t)V{0))/kBT: 



-ait/Ta 



4 Approximate Result for the VACF 

The inverse transformation of Eq.Q is only possible nu- 
merically and an exact closed form analytical expression 
for the VACF is difficult. However, an an approximate re- 
sult can be formulated using the argument of time-scale 
separation between Tc and t^. The sound waves created 
in the solvent always precedes the development of slow 



(l + 2pc/p/) 

- — ! 

a2 



Cos 



(Ol2t\ 



(16) 



with = (1 + pf/2pc) and = (1 - pj / Apl)^/'^ . For a 
neutrally buoyant particle {pc = Pf) the two contribution 
takes the form 



(v{t)vm 



2kBT 1 
3M 3^0/0 



Jo 



Ax 



g-2;t/T„^l/2 
l + a;/3 + j:2/9 



(17) 
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{vmm 



KbT _ 

e 

MI 



(3/2)(i/re) 



-\/3Si 



(18) 



To a first order, the complete VACF of a colloidal parti- 
cle can be described by a simple addition of Eq. ( [l7| and 
Eq.ra: 



3M 3 



3M 



Td 

Ivr Jo 

m 



g-3:t/r„^l/2 

l + x/3 + a;2/9 



(19) 

In Fig[5j we compare the normalized VACF of the 
Brownian particle with the theoretical predictions obtained 
by addition of Eq.(15l and Eq.([l6l). 



Finally, we compare the VACF of a Brownian parti- 
cle in the intermediate regime in Fig|6] The intermediate 
decay, between the molecular collision time and the sonic 
time, is clearly sensitive to the viscoelastic nature of the 
fluid. 




Fig. 6. Comparison of the VACF of a Brownian parti- 
cle in the intermediate decay regime. A reasonable fit to 
the data is obtained only when we consider the viscoelas- 
tic nature of the solvent. The solid line in the plot is the 
numerical evaluation of the VACF using a frequency de- 
pendent viscosity, and the dot-dashed lines is the plot of 
addition of Eq. (flSl) and Eq. (fiel) . 



In conclusion, using molecular dynamics simulations, 
we have investigated the decay of the velocity autocorre- 
lation function (VACF) of a colloid in a Lennard- Jones 
solvent. The numerical values of the shear and kinematic 
viscosities and the speed of sound, which determine the 
decay, were obtained from separate simulations of the bulk 
Lennard- Jones fluid at the same thermodynamic state point. 
These values were used to determine VACF from the ex- 
act analytical prediction. Accordingly, we divide the com- 
plete decay in three regimes, a short-time regime where 



discrete nature of the fluid plays an important role, an 
intermediate regime - governed by the interplay between 
sound propagation, vorticity diffusion and viscoelasticity 
of the fluid and a long-time regime of algebraic decay due 
to vorticity diffusion. We observe, that the decay in the 
intermediate regime can not be accurately described by 
only considering the compressibility of the fluid, but the 
viscoelastic nature of the solvent should also be taken into 
account. 
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